#include "coupling_coeff.h"
#include "Wigner_signs.h"
//#include "gsl/gsl_sf_coupling.h"
#include <fstream>
#include <iostream>
using std::cout;
using std::endl;
using std::vector;

#include <cmath>

using std::max;
using std::min;

const double INFINITE = 1E300;
// dummy variable
const int NADA = 666;

// #include <cmath>
// #include <stdlib.h>
using namespace Wigner_signs;

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ //
double CGcoeff_2(double J1_2, double m1_2, double J2_2, double m2_2, double J_2,
                 double m_2) {
  double cg = Clebsch_Gordan(J1_2 / 2.0, m1_2 / 2.0, J2_2 / 2.0, m2_2 / 2.0,
                             J_2 / 2.0, m_2 / 2.0);

  return cg;
}

double SixJSymbol_2(double J1_2, double J2_2, double J3_2, double J4_2,
                    double J5_2, double J6_2) {
  double sixJ = Wigner_6j(J1_2 / 2.0, J2_2 / 2.0, J3_2 / 2.0, J4_2 / 2.0,
                          J5_2 / 2.0, J6_2 / 2.0);

  return sixJ;
}

double NineJSymbol_2(double J1_2, double J2_2, double J3_2, double J4_2,
                     double J5_2, double J6_2, double J7_2, double J8_2,
                     double J9_2) {
  double nineJ =
      Wigner_9j(J1_2 / 2.0, J2_2 / 2.0, J3_2 / 2.0, J4_2 / 2.0, J5_2 / 2.0,
                J6_2 / 2.0, J7_2 / 2.0, J8_2 / 2.0, J9_2 / 2.0);
  return nineJ;
}
//~~~~~~~~~~~~~ start cg calculation~~~~~~~~~~~~~~~~~~~~~//

void J2M2::init(int &j2a_max) {
  this->vec.resize(j2a_max + 1);
  int size_m;
  int m2a, j2a;
  int index = 0;
  j2m2_state j2m2_t;
  for (int i = 0; i <= j2a_max; i++) {
    j2a = i;
    size_m = j2a + 1;
    this->vec[i].resize(size_m);
    for (int j = 0; j < size_m; j++) {
      m2a = 2 * j - j2a;
      this->vec[i][j] = index;
      j2m2_t.j2a = j2a;
      j2m2_t.m2a = m2a;
      j2m2_t.size = size_m;
      // j2m2_t.index = index;
      this->state.push_back(j2m2_t);
      index++;
    }
  }
}

void TrK_cg::init(J2M2 &A) {
  int size = A.state.size();
  int j2min, j2max;
  int j2a, m2a, j2b, m2b, J2, M2;
  this->vec.resize(size);
  this->exindex3.resize(size);
  this->inv_exindex3.resize(size);
  struct TrK_state state_t;
  int index = 0;
  for (int i = 0; i < size; i++) {
    this->vec[i].resize(size);
    this->exindex3[i].resize(size);
    this->inv_exindex3[i].resize(size);
    j2a = A.state[i].j2a;
    m2a = A.state[i].m2a;
    for (int j = 0; j < size; j++) {
      this->exindex3[i][j].resize(size);
      j2b = A.state[j].j2a;
      m2b = A.state[j].m2a;
      j2min = std::abs(j2a - j2b);
      j2max = j2a + j2b;
      int index_3 = 0;
      for (int k = 0; k < size; k++) {
        J2 = A.state[k].j2a;
        M2 = A.state[k].m2a;
        this->exindex3[i][j][k] = -1;
        if (M2 == (m2a + m2b) && J2 >= j2min && J2 <= j2max) {
          this->vec[i][j].push_back(index);
          // index++;
          this->exindex3[i][j][k] = index_3;
          // this->inv_exindex3[i][j].push_back(k);
          state_t.k1 = i;
          state_t.k2 = j;
          state_t.k3 = k;
          // state_t.index = index;
          index++;
          this->state.push_back(state_t);
          index_3++;
        }
      }
    }
  }
}

CG_coeff::CG_coeff() {
  /*j2a_max = 0;
  init_info = 0;*/
  int A = 0;
}
//#pragma acc routine seq
double CG_coeff::cal(int j2a, int j2b, int J2, int m2a, int m2b, int M2) const {
  // double result_3j;
  // double cg;
  int M2_x = M2 * (-1);
  // int phase;
  if ((m2a + m2b) != M2)
    return 0.0;
  if (std::abs(j2a - j2b) > J2 || (j2a + j2b) < J2)
    return 0.0;
  /*
double result_3j = gsl_sf_coupling_3j(j2a, j2b, J2, m2a, m2b, M2_x);
int phase = 1 - 2 * (int(std::abs((j2a - j2b + M2) / 2)) % 2);
double cg = result_3j * std::sqrt(J2 + 1) * phase;*/
  double cg = CGcoeff_2(j2a, m2a, j2b, m2b, J2, M2);
  return cg;
}

void CG_coeff::init_cal(J2M2 &A, TrK_cg &B) {
  // int k1, k2, k3, m;
  // int j2a, m2a, j2b, m2b, J2, M2;
  int TrK_num = B.state.size();
  // cout << TrK_num << endl;
  this->vec.resize(TrK_num);
#pragma omp parallel for
  for (int i = 0; i < TrK_num; i++) {
    int k1 = B.state[i].k1;
    int k2 = B.state[i].k2;
    int k3 = B.state[i].k3;
    int j2a = A.state[k1].j2a;
    int m2a = A.state[k1].m2a;
    int j2b = A.state[k2].j2a;
    int m2b = A.state[k2].m2a;
    int J2 = A.state[k3].j2a;
    int M2 = A.state[k3].m2a;
    this->vec[i] = this->cal(j2a, j2b, J2, m2a, m2b, M2);
  }
}
double CG_coeff::quick(int j2a, int j2b, int J2, int m2a, int m2b,
                       int M2) const {
  int m_index;
  int k1, k2, k3;
  /*
    if (J2 > j2a + j2b || J2 < abs(j2a - j2b))
      return 0.0;
    if (std::max(j2a, j2b) >= j2a_max)
      return this->cal(j2a, j2b, J2, m2a, m2b, M2);
    if (J2 >= j2a_max)
      return this->cal(j2a, j2b, J2, m2a, m2b, M2);
    m_index = (m2a + j2a) / 2;
    k1 = this->J2M2_t.vec[j2a][m_index];
    m_index = (m2b + j2b) / 2;
    k2 = this->J2M2_t.vec[j2b][m_index];

    m_index = (M2 + J2) / 2;
    k3 = this->J2M2_t.vec[J2][m_index];
    int m = this->TrKcg_t.exindex3[k1][k2][k3];
    if (m == -1) {
      // cout << " J2 > j2a+j2b or J2<abs(j2a-j2b) ? at cal_quick at
    CG_coeff::quick"
      //      << endl;
      exit(0);
      return 0.0;
    }

    int trk_num = this->TrKcg_t.vec[k1][k2][m];

    double cg = this->vec[trk_num];


    return cg; */

  return this->cal(j2a, j2b, J2, m2a, m2b, M2);
};

//~~~~~~~~~~~~~ end cg calculation~~~~~~~~~~~~~~~~~~~~~~//

//~~~~~~~~~~~~~ start 6 j calculation~~~~~~~~~~~~~~~~~~~~~~//
int JaJbJ::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

void JaJbJ::init(int &j2a_max) {
  // std::cout << "start initial JaJbJ for 6j coeff" << std::endl;
  this->vec.resize(j2a_max + 1);
  int J2max, J2min, J2num;
  int index = 0;
  for (int i = 0; i < j2a_max + 1; i++) {
    this->vec[i].resize(j2a_max + 1);
    for (int j = 0; j < j2a_max + 1; j++) {
      J2max = i + j;
      J2min = abs(i - j);
      J2num = (J2max) + 1;
      this->vec[i][j].resize(J2num);
      for (int J2 = 0; J2 < J2num; J2++)
        this->vec[i][j][J2] = -1;
      for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) { // !!!  need to check !!!
        this->vec[i][j][J2] = index;
        index++;
      }
    }
  }
  this->state.resize(index);
  index = 0;
  for (int i = 0; i <= j2a_max; i++) {
    for (int j = 0; j <= j2a_max; j++) {
      J2max = i + j;
      J2min = abs(i - j);
      for (int J2 = 0; J2 < this->vec[i][j].size(); J2++) {
        int k = this->vec[i][j][J2];
        if (k >= 0) {
          this->state[index].j2a = i;
          this->state[index].j2b = j;
          this->state[index].J2 = J2;
          this->state[index].index = k;
          index++;
        }
      }
    }
  }
}

void TrK::init(JaJbJ &A) {
  // std::cout << " start initial Trk for 6 j coeff" << std::endl;

  int size = A.state.size();
  int index = 0;
  int index_2 = 0;
  int index_3 = 0;
  TrK_state trk_state_t;

  this->vec.resize(size);
  this->state.resize(index);
  // inv_exindex2.resize(size);
  this->exindex2.resize(size);
  this->exindex3.resize(size);
  for (int i = 0; i < size; i++) {
    exindex2[i].resize(size);
    index_2 = 0;
    for (int j = 0; j < size; j++) {
      this->exindex2[i][j] = -1;
      if (A.state[i].j2b == A.state[j].j2a) {
        this->exindex2[i][j] = index_2;
        index_2++;
      }
    }
    int size2 = index_2;
    this->exindex3[i].resize(size2);
    this->vec[i].resize(size2);
    int index_22 = -1;
    for (int j = 0; j < size; j++) {
      if (A.state[i].j2b == A.state[j].j2a) {
        index_22++;
        this->exindex3[i][index_22].resize(size);
        index_3 = 0;
        for (int k = 0; k < size; k++) {
          this->exindex3[i][index_22][k] = -1;
          if (A.state[k].j2a == A.state[i].J2 &&
              A.state[k].j2b == A.state[j].j2b) {
            exindex3[i][index_22][k] = index_3;
            this->vec[i][index_22].push_back(index);
            trk_state_t.k1 = i;
            trk_state_t.k2 = j;
            trk_state_t.k3 = k;
            // trk_state_t.index = index;
            this->state.push_back(trk_state_t);
            index_3++;
            index++;
          }
        }
      }
    }
  }
  // cout << sizeof inv_exindex3 << endl;
  /*
  this->vec.resize(size);
  for (int i = 0; i < size; i++) {
    int size_2 = inv_exindex2[i].size();
    this->vec[i].resize(size_2);
    for (int j = 0; j < size_2; j++) {
      int size_3 = inv_exindex3[i][j].size();
      this->vec[i][j].resize(size_3);
      for (int k = 0; k < size_3; k++) {
        this->vec[i][j][k] = index;
        index++;
      }
    }
  }
  this->state.resize(index);
  index = 0;
  for (int i = 0; i < size; i++) {
    int size_2 = inv_exindex2[i].size();
    for (int j = 0; j < size_2; j++) {
      int size_3 = inv_exindex3[i][j].size();
      for (int k = 0; k < size_3; k++) {
        this->state[index].k1 = i;
        this->state[index].k2 = inv_exindex2[i][j];
        this->state[index].k3 = inv_exindex3[i][j][k];
        this->state[index].index = index;
        index++;
      }
    }
  }*/
}

SixJ_coeff::SixJ_coeff() {
  // j2a_max = 0;
  // init_info = -1;
  int n = 0;
  // JaJbJ_t.init(n);
}

double SixJ_coeff::cal(int j2a, int j2b, int J2ab, int j2c, int J2,
                       int J2bc) const {
  return SixJSymbol_2(j2a, j2b, J2ab, j2c, J2, J2bc);
  // return gsl_sf_coupling_6j(j2a, j2b, J2ab, j2c, J2, J2bc);
}

void SixJ_coeff::init_cal(JaJbJ &A, TrK &B) {
  // std::cout << " start initial 6 j cal" << std::endl;
  int size = B.state.size();
  // int j2a, j2b, J2ab, j2c, J2, J2bc;
  // int j2b_t, J2ab_t, j2c_t;
  this->vec.resize(size);
#pragma omp parallel for
  for (int i = 0; i < size; i++) {
    int k1 = B.state[i].k1;
    int k2 = B.state[i].k2;
    int k3 = B.state[i].k3;
    int j2a = A.state[k1].j2a;
    int j2b = A.state[k1].j2b;
    int J2ab = A.state[k1].J2;
    int j2b_t = A.state[k2].j2a;
    int j2c = A.state[k2].j2b;
    int J2bc = A.state[k2].J2;
    int J2ab_t = A.state[k3].j2a;
    int j2c_t = A.state[k3].j2b;
    int J2 = A.state[k3].J2;
    if (j2b_t != j2b || j2c_t != j2c || J2ab_t != J2ab) {
      std::cerr << "error at init_cal of sixJ_coeff" << endl;
      exit(0);
    }
    // double sixj_gsl = gsl_sf_coupling_6j(j2a, j2b, J2ab, j2c, J2, J2bc);
    // cout << j2a << "\t" << j2b << "\t" << J2ab << "\t" << j2c << "\t" << J2
    //      << "\t" << J2bc << "\t"<<sixj_t<<"\t "<<endl;

    double sixj = this->cal(j2a, j2b, J2ab, j2c, J2, J2bc);
    this->vec[i] = sixj;
    // if (std::abs(sixj_gsl - sixj) > 0.0001) {
    //   cout << j2a << "\t" << j2b << "\t" << J2ab << "\t" << j2c << "\t" << J2
    //        << "\t" << J2bc << "\t sixj_gsl = " << sixj_gsl
    //        << "\t sixj = " << sixj << endl;
    //   exit(0);
    // }
    // cout <<"\t -- "<< sixj << endl;
  }
  cout << "Finished Init SixJ" << endl;
}

double SixJ_coeff::quick(int j2a, int j2b, int J2ab, int j2c, int J2,
                         int J2bc) const {
  /*int tri1 = this->JaJbJ_t.Tri(j2a, j2b, J2ab);
  int tri2 = this->JaJbJ_t.Tri(j2b, j2c, J2bc);
  int tri3 = this->JaJbJ_t.Tri(J2ab, j2c, J2);
  int tri4 = this->JaJbJ_t.Tri(j2a, J2bc, J2);
  */
  // if (tri1 != 1 || tri2 != 1 || tri3 != 1 || tri4 != 1) {
  // return 0.0;
  //}
  /*
  if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
    return 0.0;
  if (J2bc < std::abs(j2b - j2c) || J2bc > (j2b + j2c))
    return 0.0;
  if (J2 < std::abs(J2ab - j2c) || J2 > (J2ab + j2c))
    return 0.0;
  if (std::max(j2a, j2b) >= this->j2a_max ||
      std::max(J2ab, j2c) >= this->j2a_max || J2bc >= this->j2a_max) {
    // return gsl_sf_coupling_6j(j2a, j2b, J2ab, j2c, J2, J2bc);
    return this->cal(j2a, j2b, J2ab, j2c, J2, J2bc);
    // exit(0);
  }

  if (std::max(j2a, j2b) > this->j2a_max ||
      std::max(j2b, j2c) > this->j2a_max) {
    // std::cerr << "wrong at sixJ_quick need to set a larger j2a_max" <<
    // std::endl;
    exit(0);
  } else {
    int k1 = this->JaJbJ_t.vec[j2a][j2b][J2ab];
    int k2 = this->JaJbJ_t.vec[j2b][j2c][J2bc];
    int k3 = this->JaJbJ_t.vec[J2ab][j2c][J2];
    if (k1 < 0 || k2 < 0 || k3 < 0) {
      // cout << "wrong at sixJ_quick k <0 " << endl;
      // cout << "k"
      //      << "\t j2a"
      //      << "\t j2b"
      //      << "\t J2ab" << endl;
      // cout << k1 << "\t" << j2a << "\t" << j2b << "\t" << J2ab << endl;
      // cout << "k"
      //      << "\t j2b"
      //      << "\t j2c"
      //      << "\t J2bc" << endl;
      // cout << k2 << "\t" << j2b << "\t" << j2c << "\t" << J2bc << endl;
      // cout << "k"
      //      << "\t J2ab"
      //      << "\t j2c"
      //      << "\t J2" << endl;
      // cout << k3 << "\t" << J2ab << "\t" << j2c << "\t" << J2 << endl;
      exit(0);
    }
    int m = this->TrK_t.exindex2[k1][k2];
    int n = this->TrK_t.exindex3[k1][m][k3];
    int trk = this->TrK_t.vec[k1][m][n];
    double sixj = this->vec[trk];
    return sixj;
  }*/
  return this->cal(j2a, j2b, J2ab, j2c, J2, J2bc);
}
//~~~~~~~~~~~~~ end 6 j calculation~~~~~~~~~~~~~~~~~~~~~~//

//~~~~~~~~~~~~ start 9j calculation~~~~~~~~~~~~~~~~~~~~~~//

void TrK_ninJ::init(JaJbJ &A) {
  // cout << "start to init TrK_ninJ for 9j" << endl;
  int size = A.state.size();
  this->vec.resize(size);
  this->exindex3.resize(size);
  TrK_state state_t;
  int index = 0;
  for (int i = 0; i < size; i++) {
    int j2a = A.state[i].j2a;
    int j2b = A.state[i].j2b;
    int J2ab = A.state[i].J2;
    this->vec[i].resize(size);
    this->exindex3[i].resize(size);
    for (int j = 0; j < size; j++) {
      int j2c = A.state[j].j2a;
      int j2d = A.state[j].j2b;
      int J2cd = A.state[j].J2;
      this->exindex3[i][j].resize(size);
      int index_3 = 0;
      for (int k = 0; k < size; k++) {
        int J2ac = A.state[k].j2a;
        int J2bd = A.state[k].j2b;
        int J2 = A.state[k].J2;
        int tri1 = A.Tri(j2a, j2c, J2ac);
        int tri2 = A.Tri(j2b, j2d, J2bd);
        int tri3 = A.Tri(J2ab, J2cd, J2);
        this->exindex3[i][j][k] = -1;
        if (tri1 > 0 && tri2 > 0 && tri3 > 0) {
          this->vec[i][j].push_back(index);
          this->exindex3[i][j][k] = index_3;
          state_t.k1 = i;
          state_t.k2 = j;
          state_t.k3 = k;
          this->state.push_back(state_t);
          index++;
          index_3++;
        }
      }
    }
  }
}
NinJ_coeff::NinJ_coeff() {
  /*j2a_max = 0;
  init_info = 0;
  SixJ_coeff sixJ_t;
  sixJ_t.init(0);
  sixJ = sixJ_t;*/
  int A = 0;
}

double NinJ_coeff::cal(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
                       int J2ac, int J2bd, int J2) const {
  return NineJSymbol_2(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
  // return gsl_sf_coupling_9j(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
}

void NinJ_coeff::init_cal(JaJbJ &A, TrK_ninJ &B) {
  // cout << "start to initial init_cal for 9j " << endl;
  int size = B.state.size();
  this->vec.resize(size);
  // cout << " 9j calculation times " << size << endl;
#pragma omp parallel for
  for (int i = 0; i < size; i++) {
    int k1 = B.state[i].k1;
    int k2 = B.state[i].k2;
    int k3 = B.state[i].k3;
    int j2a = A.state[k1].j2a;
    int j2b = A.state[k1].j2b;
    int J2ab = A.state[k1].J2;
    int j2c = A.state[k2].j2a;
    int j2d = A.state[k2].j2b;
    int J2cd = A.state[k2].J2;
    int J2ac = A.state[k3].j2a;
    int J2bd = A.state[k3].j2b;
    int J2 = A.state[k3].J2;
    // this->vec[i] =
    //     gsl_sf_coupling_9j(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
    this->vec[i] = this->cal(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
  }
}

// double NinJ_coeff::quick2(int j2a, int j2b, int J2ab, int j2c, int j2d,
//                           int J2cd, int J2ac, int J2bd, int J2) const {
//   /*int tri1 = this->JaJbJ_t.Tri(j2a, j2b, J2ab);
//   int tri2 = this->JaJbJ_t.Tri(j2c, j2d, J2cd);
//   int tri3 = this->JaJbJ_t.Tri(J2ab, J2cd, J2);
//   int tri4 = this->JaJbJ_t.Tri(J2ac, J2bd, J2);
//   int tri5 = this->JaJbJ_t.Tri(j2a, j2c, J2ac);
//   int tri6 = this->JaJbJ_t.Tri(j2b, j2d, J2bd);
//   if (tri1 < 0 || tri2 < 0 || tri3 < 0)
//     return 0.0;
//   if (tri4 < 0 || tri5 < 0 || tri6 < 0)
//     return 0.0;
//     */
//   if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
//     return 0.0;
//   if (J2cd < std::abs(j2c - j2d) || J2cd > (j2c + j2d))
//     return 0.0;
//   if (J2 < std::abs(J2ab - J2cd) || J2 > (J2ab + J2cd))
//     return 0.0;
//   if (J2ac < std::abs(j2a - j2c) || J2ac > (j2a + j2c))
//     return 0.0;
//   if (J2bd < std::abs(j2b - j2d) || J2bd > (j2b + j2d))
//     return 0.0;
//   if (J2 < std::abs(J2ac - J2bd) || J2 > (J2ac + J2bd))
//     return 0.0;
//   if (std::max(j2a, j2b) >= this->j2a_max ||
//       std::max(j2c, j2d) >= this->j2a_max)
//     // return gsl_sf_coupling_9j(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd,
//     // J2);
//     return this->cal(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
//   if (std::max(J2ab, J2cd) >= this->j2a_max ||
//       std::max(J2ac, J2bd) >= this->j2a_max || J2 >= this->j2a_max)
//     // return gsl_sf_coupling_9j(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd,
//     // J2);
//     return this->cal(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
//   int k1 = this->JaJbJ_t.vec[j2a][j2b][J2ab];
//   int k2 = this->JaJbJ_t.vec[j2c][j2d][J2cd];
//   int k3 = this->JaJbJ_t.vec[J2ac][J2bd][J2];
//   int m = this->TrKninJ_t.exindex3[k1][k2][k3];
//   if (m < 0) {
//     cout << "9j condit disapprove at ninJ_coeff :: quick  m < 0" << endl;
//     return 0.0;
//   }
//   // cout << k1 << "\t" << k2 << "\t" << k3 << "\t" << m << endl;
//   int tri_num = this->TrKninJ_t.vec[k1][k2][m];
//   double ninJ = this->vec[tri_num];
//   return ninJ;
// }
double NinJ_coeff::quick(int &j2a, int &j2b, int &J2ab, int &j2c, int &j2d,
                         int &J2cd, int &J2ac, int &J2bd, int &J2) const {

  // http://mathworld.wolfram.com/Wigner9j-Symbol.html (3)
  return NineJSymbol_2(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
  //
  // if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
  //   return 0.0;
  // if (J2cd < std::abs(j2c - j2d) || J2cd > (j2c + j2d))
  //   return 0.0;
  // if (J2 < std::abs(J2ab - J2cd) || J2 > (J2ab + J2cd))
  //   return 0.0;
  // if (J2ac < std::abs(j2a - j2c) || J2ac > (j2a + j2c))
  //   return 0.0;
  // if (J2bd < std::abs(j2b - j2d) || J2bd > (j2b + j2d))
  //   return 0.0;
  // if (J2 < std::abs(J2ac - J2bd) || J2 > (J2ac + J2bd))
  //   return 0.0;
  // int g2_min = std::max(std::abs(j2b - J2cd), std::abs(j2c - J2bd));
  // g2_min = std::max(g2_min, std::abs(J2 - j2a));
  // int g2_max = std::min(j2b + J2cd, j2c + J2bd);
  // g2_max = std::min(g2_max, J2 + j2a);
  //
  // double ninJ = 0.0;
  // for (int g2 = g2_min; g2 <= g2_max; g2 = g2 + 2) {
  //   double part1_t = 0.0;
  //   part1_t = sixJ.quick(j2a, j2b, J2ab, J2cd, J2, g2);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * sixJ.quick(j2c, j2d, J2cd, j2b, g2, J2bd);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * sixJ.quick(J2ac, J2bd, J2, g2, j2a, j2c);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * phase(g2) * (g2 + 1);
  //   ninJ += part1_t;
  // }
  // return ninJ;
}
double NinJ_coeff::quick(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
                         int J2ac, int J2bd, int J2) const {

  // http://mathworld.wolfram.com/Wigner9j-Symbol.html (3)
  return NineJSymbol_2(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd, J2);
  //
  // if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
  //   return 0.0;
  // if (J2cd < std::abs(j2c - j2d) || J2cd > (j2c + j2d))
  //   return 0.0;
  // if (J2 < std::abs(J2ab - J2cd) || J2 > (J2ab + J2cd))
  //   return 0.0;
  // if (J2ac < std::abs(j2a - j2c) || J2ac > (j2a + j2c))
  //   return 0.0;
  // if (J2bd < std::abs(j2b - j2d) || J2bd > (j2b + j2d))
  //   return 0.0;
  // if (J2 < std::abs(J2ac - J2bd) || J2 > (J2ac + J2bd))
  //   return 0.0;
  // int g2_min = std::max(std::abs(j2b - J2cd), std::abs(j2c - J2bd));
  // g2_min = std::max(g2_min, std::abs(J2 - j2a));
  // int g2_max = std::min(j2b + J2cd, j2c + J2bd);
  // g2_max = std::min(g2_max, J2 + j2a);
  //
  // double ninJ = 0.0;
  // for (int g2 = g2_min; g2 <= g2_max; g2 = g2 + 2) {
  //   double part1_t = 0.0;
  //   part1_t = sixJ.quick(j2a, j2b, J2ab, J2cd, J2, g2);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * sixJ.quick(j2c, j2d, J2cd, j2b, g2, J2bd);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * sixJ.quick(J2ac, J2bd, J2, g2, j2a, j2c);
  //   if (part1_t == 0.0)
  //     continue;
  //   part1_t = part1_t * phase(g2) * (g2 + 1);
  //   ninJ += part1_t;
  // }
  // return ninJ;
}

//~~~~~~~~~~~~ end 9j calculation~~~~~~~~~~~~~~~~~~~~~~//

//~~~~~~~~~ start  binomial coeff ~~~~~~~~~~~//

double Binom::fac_cal(int n) {
  // const int n = make_int(x);
  //#pragma acc data
  // copyin(fact_table_coupling[0:170],doub_fact_table_coupling[0:297])
  double fact = 1.0; // this->fac2vec[n];
  // int gsl_fac_max = 170;
  // if (n < gsl_fac_max) {
  //   fact = fact_table_coupling[n].f;
  //   return fact;
  // } else {
  for (int i = 1; i <= n; i++)
    fact *= i;

  return fact;
  // }
}
double Binom::fac2_cal(int n) {
  double double_fact = 1.0; // this->fac2vec[n];
  // int gsl_fac2_max = 297;
  // if (n < gsl_fac2_max) {
  //   double_fact = doub_fact_table_coupling[n].f;
  //   return double_fact;
  // } else {
  for (int i = n; i > 0; i -= 2)
    double_fact *= i;
  return double_fact;
  // }
}
void Binom::facveccal(int &n) {
  // n*(n-2)*(n-4)...
  this->facvec.resize(n + 1);
  this->facvec[0] = 1;
  this->facvec[1] = 1;
  for (int i = 2; i < n; i++) {
    double k = i;
    for (int j = i - 1; j > 0; j = j - 1) {
      k = k * j;
    }
    this->facvec[i] = k;
  }
}
void Binom::fac2veccal(int &n) {
  // n*(n-2)*(n-4)...
  this->fac2vec.resize(n + 2);
  this->fac2vec[0] = 1;
  this->fac2vec[1] = 1;
  for (int i = 2; i < n; i++) {
    double k = i;
    for (int j = i - 2; j > 0; j = j - 2) {
      k = k * j;
    }
    this->fac2vec[i] = k;
  }
}
double Binom::binom_cal(int n, int k1, int k2) {
  // if (n < 0 || k1 > n || k2 > n) {
  //   //cout << "k > n in binomial binom" << endl;
  //   exit(0);
  // }
  double C_n = this->fac2_cal(n);
  double C_k1 = this->fac2_cal(k1);
  double C_k2 = this->fac2_cal(k2);
  // double bino =
  //     this->fac2vec[n] / (1.0 * this->fac2vec[k1] * this->fac2vec[k2]);
  double temp = 1.0 * C_k1 * C_k2;
  double bino = C_n / temp;
  return bino;
}

double Binom::init_cal(int &n_max) {
  this->vec.resize(n_max + 1);
  for (int i = 0; i < n_max; i++) {
    this->vec[i].resize(n_max + 1);
    int size2 = std::floor(i + 1);
    for (int j = 0; j < size2; j++) {
      this->vec[i][j].resize(n_max + 1);
      for (int k = 0; k < size2; k++) {
        this->vec[i][j][k] = this->binom_cal(i, j, k);
        // this->vec[i][i - j] = this->binom_cal(i, j);
      }
    }
  }
  return 0.0;
}

//~~~~~~~~~ end  binomial coeff ~~~~~~~~~~~//

//~~~~~~~~~ begin  harmonic-oscillator brackets ~~~~~~~~~~~//
// Hobra::Hobra(){
//   Binom binom_t;
//   int n = 0;
//   binom_t.init(n);
//   binom_p =binom_t;
//   NinJ_coeff ninJ_t;
//   ninJ_p = ninJ_t;
//   CG_coeff cg_t;
//   cg_p = cg_t;
//
// }

void Hobra::facveccal(int n) {
  // n*(n-2)*(n-4)...
  this->facvec.resize(n + 2);
  this->facvec[0] = 1;
  this->facvec[1] = 1;
  for (int i = 2; i < n; i++) {
    double k = i;
    for (int j = i - 1; j > 0; j = j - 1) {
      k = k * j;
    }
    this->facvec[i] = k;
  }
}
void Hobra::fac2veccal(int n) {
  // n*(n-2)*(n-4)...
  this->fac2vec.resize(n + 2);
  this->fac2vec[0] = 1;
  this->fac2vec[1] = 1;
  for (int i = 2; i < n; i++) {
    double k = i;
    for (int j = i - 2; j > 0; j = j - 2) {
      k = k * j;
    }
    this->fac2vec[i] = k;
  }
}
double Hobra::binom_cal(int n, int k1, int k2) const {
  // if (n < 0 || k1 > n || k2 > n) {
  //   //cout << "k > n in binomial binom" << endl;
  //   exit(0);
  // }
  double bino =
      this->fac2vec[n] / (1.0 * this->fac2vec[k1] * this->fac2vec[k2]);
  return bino;
}
void Hobra::init_binom(int n_max_fac, int n_max_fac_2) {
  this->facveccal(n_max_fac);
  this->fac2veccal(n_max_fac_2);
  this->vec.resize(n_max_fac + 1);
  for (int i = 0; i < n_max_fac; i++) {
    this->vec[i].resize(n_max_fac_2 + 1);
    int size2 = std::floor(i + 1);
    for (int j = 0; j < n_max_fac_2; j++) {
      this->vec[i][j].resize(n_max_fac_2 + 1);
      for (int k = 0; k < n_max_fac_2; k++) {
        this->vec[i][j][k] = this->binom_cal(i, j, k);
        // this->vec[i][i - j] = this->binom_cal(i, j);
      }
    }
  }
}

double Hobra::fac_cal(int n) {
  // const int n = make_int(x);
  double fact = 1.0;
  if (n < 160) {
    fact = this->facvec[n];
    return fact;
  } else {

    for (int i = 1; i <= n; i++)
      fact *= i;

    return fact;
  }
}
double Hobra::fac2_cal(int n) {
  double double_fact = 1.0;
  if (n < 180) {
    double_fact = this->fac2vec[n];
    return double_fact;
  } else {
    for (int i = n; i > 0; i -= 2)
      double_fact *= i;
    return double_fact;
  }
}

double Hobra::G(int &e1, int &l1, int &e_a, int &l_a, int &e_b,
                int &l_b) const {
  // double value_g;
  double part1, part2, part3;
  int l2a = l_a * 2, l2b = l_b * 2, l21 = l1 * 2;
  part1 = std::sqrt((l2a + 1) * (l2b + 1));
  // double result_3j = gsl_sf_coupling_3j(l2a, l2b, l21, 0, 0, 0);
  // int phase = 1 - 2 * (int(std::abs((l2a - l2b) / 2)) % 2);
  // double cg = result_3j * std::sqrt(l21 + 1) * phase;
  // part2 = cg;
  // return 0.0;
  part2 = this->cg_p.quick(l2a, l2b, l21, 0, 0, 0);
  if (part2 == 0)
    return 0.0;
  int n2_1, n2_2, n2_3;
  n2_1 = e1 - l1;
  n2_2 = e_a - l_a;
  n2_3 = e_b - l_b;
  // if (n2_1 < 0 || n2_2 < 0 || n2_3 < 0) {
  //   //std::cout << "wrong at HObra G , !! 2n = e - l < 0" << endl;
  //   exit(0);
  // }
  double temp_1 = binom_p.binom_cal(n2_1, n2_2, n2_3);
  double temp_2 = binom_p.binom_cal(e1 + l1 + 1, e_a + l_a + 1, e_b + l_b + 1);
  part3 = temp_1 * temp_2;
  // part3 = this->binom_p.vec[n2_1][n2_2][n2_3] *
  //         this->binom_p.vec[e1 + l1 + 1][e_a + l_a + 1][e_b + l_b + 1];
  double value_g = part1 * part2 * std::sqrt(part3);
  return value_g;
  // return value_g;
}

double Hobra::G_2(int &e1, int &l1, int &e_a, int &l_a, int &e_b,
                  int &l_b) const {
  // double value_g;
  // return 1.0;
  double part1, part2, part3;
  int l2a = l_a * 2, l2b = l_b * 2, l21 = l1 * 2;
  part1 = std::sqrt((l2a + 1) * (l2b + 1));
  // double result_3j = gsl_sf_coupling_3j(l2a, l2b, l21, 0, 0, 0);
  // int phase = 1 - 2 * (int(std::abs((l2a - l2b) / 2)) % 2);
  // double cg = result_3j * std::sqrt(l21 + 1) * phase;
  // part2 = cg;
  // return 0.0;
  part2 = this->cg_p.quick(l2a, l2b, l21, 0, 0, 0);
  if (part2 == 0)
    return 0.0;
  int n2_1, n2_2, n2_3;
  n2_1 = e1 - l1;
  n2_2 = e_a - l_a;
  n2_3 = e_b - l_b;
  // if (n2_1 < 0 || n2_2 < 0 || n2_3 < 0) {
  //   //std::cout << "wrong at HObra G , !! 2n = e - l < 0" << endl;
  //   exit(0);
  // }
  double temp_1 = this->binom_cal(n2_1, n2_2, n2_3);
  double temp_2 = this->binom_cal(e1 + l1 + 1, e_a + l_a + 1, e_b + l_b + 1);
  part3 = temp_1 * temp_2;
  // part3 = this->binom_p.vec[n2_1][n2_2][n2_3] *
  //         this->binom_p.vec[e1 + l1 + 1][e_a + l_a + 1][e_b + l_b + 1];
  double value_g = part1 * part2 * std::sqrt(part3);
  return value_g;
  // return value_g;
}

double Hobra::cal(int &E, int &L, int &e, int &l, int &e1, int &l1, int &e2,
                  int &l2, int &lambda, double &d) const {
  double part1;
  double t1, t2;
  if ((E + e) != (e1 + e2))
    return 0.0;
  if (lambda > (l + L) || lambda < std::abs(l - L))
    return 0.0;
  if (lambda > (l1 + l2) || lambda < std::abs(l1 - l2))
    return 0.0;

  int N = (E - L) / 2;
  int n = (e - l) / 2;
  int n1 = (e1 - l1) / 2;
  int n2 = (e2 - l2) / 2;
  int phase_t =
      N + n + n1 + n2; // from the different convension of H.O. Wave_Fun
  int phase_diff = phase_cal(phase_t);

  int e_a, e_b, e_c, e_d, l_a, l_b, l_c, l_d;
  t1 = std::pow(d, (e1 - e));
  t2 = std::pow((1 + d), (e1 + e2));
  part1 = std::sqrt(t1 / (1.0 * t2));
  int e_d_max = std::min(e, e2);
  double part2 = 0.0;
  for (e_d = 0; e_d <= e_d_max; e_d++) {
    e_a = e1 - e + e_d;
    e_b = e - e_d;
    e_c = e2 - e_d;
    double part2_1 = std::pow(-1 * d, e_d);
    double part2_2 = 0;
    for (l_a = e_a; l_a >= 0; l_a = l_a - 2) {
      for (l_b = e_b; l_b >= 0; l_b = l_b - 2) {
        if (l1 > (l_a + l_b) || l1 < std::abs(l_a - l_b))
          continue;
        double G1 = this->G(e1, l1, e_a, l_a, e_b, l_b);
        for (l_c = e_c; l_c >= 0; l_c = l_c - 2) {
          if (L > (l_a + l_c) || L < std::abs(l_a - l_c))
            continue;
          double G3 = this->G(E, L, e_a, l_a, e_c, l_c);
          for (l_d = e_d; l_d >= 0; l_d = l_d - 2) {
            if (l > (l_b + l_d) || l < std::abs(l_b - l_d))
              continue;
            if (l2 > (l_c + l_d) || l2 < std::abs(l_c - l_d))
              continue;

            double part2_t =
                this->ninJ_p.quick(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d * 2,
                                   l2 * 2, L * 2, l * 2, lambda * 2);
            // double part2_t =
            //     gsl_sf_coupling_9j(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d *
            //     2,
            //                        l2 * 2, L * 2, l * 2, lambda * 2);
            if (part2_t != 0) {
              // double G2 = this->G(e2, l2, e_c, l_c, e_d, l_d);
              // double G4 = this->G(e, l, e_b, l_b, e_d, l_d);
              part2_t = part2_t * G1 * G3 *
                        this->G(e2, l2, e_c, l_c, e_d, l_d) *
                        this->G(e, l, e_b, l_b, e_d, l_d);
            }
            part2_2 += part2_t;
          }
        }
      }
    }
    part2 += part2_1 * part2_2;
  }
  double value = part1 * part2 * phase_diff;
  return value;
}

double Hobra::cal2(int &E, int &L, int &e, int &l, int &e1, int &l1, int &e2,
                   int &l2, int &lambda, double &d) const {
  double part1;
  double t1, t2;
  if ((E + e) != (e1 + e2))
    return 0.0;
  if (lambda > (l + L) || lambda < std::abs(l - L))
    return 0.0;
  if (lambda > (l1 + l2) || lambda < std::abs(l1 - l2))
    return 0.0;

  int N = (E - L) / 2;
  int n = (e - l) / 2;
  int n1 = (e1 - l1) / 2;
  int n2 = (e2 - l2) / 2;
  int phase_t =
      N + n + n1 + n2; // from the different convension of H.O. Wave_Fun
  int phase_diff = phase_cal(phase_t);

  int e_a, e_b, e_c, e_d, l_a, l_b, l_c, l_d;
  t1 = std::pow(d, (e1 - e));
  t2 = std::pow((1 + d), (e1 + e2));
  part1 = std::sqrt(t1 / (1.0 * t2));
  int e_d_max = std::min(e, e2);
  double part2 = 0.0;
  for (e_d = 0; e_d <= e_d_max; e_d++) {
    e_a = e1 - e + e_d;
    e_b = e - e_d;
    e_c = e2 - e_d;
    double part2_1 = std::pow(-1 * d, e_d);
    double part2_2 = 0;
    for (l_a = e_a; l_a >= 0; l_a = l_a - 2) {
      for (l_b = e_b; l_b >= 0; l_b = l_b - 2) {
        if (l1 > (l_a + l_b) || l1 < std::abs(l_a - l_b))
          continue;
        double G1 = this->G_2(e1, l1, e_a, l_a, e_b, l_b);
        for (l_c = e_c; l_c >= 0; l_c = l_c - 2) {
          if (L > (l_a + l_c) || L < std::abs(l_a - l_c))
            continue;
          double G3 = this->G_2(E, L, e_a, l_a, e_c, l_c);
          for (l_d = e_d; l_d >= 0; l_d = l_d - 2) {
            if (l > (l_b + l_d) || l < std::abs(l_b - l_d))
              continue;
            if (l2 > (l_c + l_d) || l2 < std::abs(l_c - l_d))
              continue;

            double part2_t =
                this->ninJ_p.quick(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d * 2,
                                   l2 * 2, L * 2, l * 2, lambda * 2);
            // double part2_t =
            //     gsl_sf_coupling_9j(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d *
            //     2,
            //                        l2 * 2, L * 2, l * 2, lambda * 2);
            if (part2_t != 0) {
              // double G2 = this->G(e2, l2, e_c, l_c, e_d, l_d);
              // double G4 = this->G(e, l, e_b, l_b, e_d, l_d);
              part2_t = part2_t * G1 * G3 *
                        this->G_2(e2, l2, e_c, l_c, e_d, l_d) *
                        this->G_2(e, l, e_b, l_b, e_d, l_d);
            }
            part2_2 += part2_t;
          }
        }
      }
    }
    part2 += part2_1 * part2_2;
  }
  double value = part1 * part2 * phase_diff;
  return value;
}

//~~~~~~~~~ end  harmonic-oscillator brackets ~~~~~~~~~~~//
// ~~~~~~~~~  Pre_coeff

double Pre_coeff::fac_cal_this(int n) {
  double fact = 1.0; // this->fac2vec[n];
  for (int i = 1; i <= n; i++)
    fact *= i;

  return fact;
}

double Pre_coeff::fac2_cal_this(int n) {
  double double_fact = 1.0; // this->fac2vec[n];
  for (int i = n; i > 0; i -= 2)
    double_fact *= i;
  return double_fact;
}

double Pre_coeff::fac_cal(int n) {
  // const int n = make_int(x);
  double fact = 1.0;
  if (n < 120) {
    fact = this->facvec[n];
    return fact;
  } else {

    for (int i = 1; i <= n; i++)
      fact *= i;

    return fact;
  }
}
double Pre_coeff::fac2_cal(int n) {
  double double_fact = 1.0;
  if (n < 120) {
    double_fact = this->fac2vec[n];
    return double_fact;
  } else {
    for (int i = n; i > 0; i -= 2)
      double_fact *= i;
    return double_fact;
  }
}

double Pre_coeff::binom_cal_this(int n, int k1, int k2) {
  // return 1.0;//
  double C_n = this->fac2_cal_this(n);
  double C_k1 = this->fac2_cal_this(k1);
  double C_k2 = this->fac2_cal_this(k2);
  // double bino =
  //     this->fac2vec[n] / (1.0 * this->fac2vec[k1] * this->fac2vec[k2]);
  double temp = 1.0 * C_k1 * C_k2;
  double bino = C_n / temp;
  return bino;
}

void Pre_coeff::init_Binom() {
  int n_max = 150;
  this->facvec.resize(n_max + 1);
  this->facvec[0] = 1;
  this->facvec[1] = 1;

  this->fac2vec.resize(n_max + 1);
  this->fac2vec[0] = 1;
  this->fac2vec[1] = 1;
  for (int i = 2; i < n_max; i++) {
    this->facvec[i] = fac_cal_this(i);
    this->fac2vec[i] = fac2_cal_this(i);
  }

  this->vec_bino.resize(n_max + 1);
  for (int i = 0; i < n_max; i++) {
    this->vec_bino[i].resize(n_max + 1);
    int size2 = std::floor(i + 1);
    for (int j = 0; j < n_max; j++) {
      this->vec_bino[i][j].resize(n_max + 1);
      for (int k = 0; k < n_max; k++) {
        this->vec_bino[i][j][k] = this->binom_cal_this(i, j, k);
        // this->vec[i][i - j] = this->binom_cal(i, j);
      }
    }
  }
}

double Pre_coeff::cg_cal_this(int j2a, int j2b, int J2, int m2a, int m2b,
                              int M2) {
  // int M2_x = M2 * (-1);
  // int phase;
  if ((m2a + m2b) != M2)
    return 0.0;
  if (std::abs(j2a - j2b) > J2 || (j2a + j2b) < J2)
    return 0.0;
  double cg = CGcoeff_2(j2a, m2a, j2b, m2b, J2, M2);
  return cg;
}

double Pre_coeff::cg_quick(int j2a, int j2b, int J2, int m2a, int m2b, int M2) {
  int m_index;
  int k1, k2, k3;
  return cg_cal_this(j2a, j2b, J2, m2a, m2b, M2);
  if (J2 > j2a + j2b || J2 < abs(j2a - j2b))
    return 0.0;
  if (std::max(j2a, j2b) >= j2a_max_cg)
    return cg_cal_this(j2a, j2b, J2, m2a, m2b, M2);
  if (J2 >= j2a_max_cg)
    return cg_cal_this(j2a, j2b, J2, m2a, m2b, M2);
  m_index = (m2a + j2a) / 2;
  k1 = this->J2M2_t.vec[j2a][m_index];
  m_index = (m2b + j2b) / 2;
  k2 = this->J2M2_t.vec[j2b][m_index];

  m_index = (M2 + J2) / 2;
  k3 = this->J2M2_t.vec[J2][m_index];
  int m = this->TrKcg_t.exindex3[k1][k2][k3];
  if (m == -1) {
    // cout << " J2 > j2a+j2b or J2<abs(j2a-j2b) ? at cal_quick at
    // CG_coeff::quick"
    //      << endl;
    // exit(0);
    return 0.0;
  }

  int trk_num = this->TrKcg_t.vec[k1][k2][m];

  double cg = this->vec_cg[trk_num];
  return cg;
}
//
double Pre_coeff::sixJ_this(int j2a, int j2b, int J2ab, int j2c, int J2,
                            int J2bc) {
  return SixJSymbol_2(j2a, j2b, J2ab, j2c, J2, J2bc);
}

double Pre_coeff::SixJ_quick(int j2a, int j2b, int J2ab, int j2c, int J2,
                             int J2bc) {
  double sixj = 0.0; //
  int tri1 = this->JaJbJ_t.Tri(j2a, j2b, J2ab);
  int tri2 = this->JaJbJ_t.Tri(j2b, j2c, J2bc);
  int tri3 = this->JaJbJ_t.Tri(J2ab, j2c, J2);
  int tri4 = this->JaJbJ_t.Tri(j2a, J2bc, J2);

  if (tri1 != 1 || tri2 != 1 || tri3 != 1 || tri4 != 1) {
    return 0.0;
  }

  if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
    return 0.0;
  if (J2bc < std::abs(j2b - j2c) || J2bc > (j2b + j2c))
    return 0.0;
  if (J2 < std::abs(J2ab - j2c) || J2 > (J2ab + j2c))
    return 0.0;
  if (std::max(j2a, j2b) >= this->j2a_max_6j ||
      std::max(J2ab, j2c) >= this->j2a_max_6j || J2bc >= this->j2a_max_6j) {
    // return gsl_sf_coupling_6j(j2a, j2b, J2ab, j2c, J2, J2bc);
    return sixJ_this(j2a, j2b, J2ab, j2c, J2, J2bc);
    // exit(0);
  }

  if (std::max(j2a, j2b) > this->j2a_max_6j ||
      std::max(j2b, j2c) > this->j2a_max_6j) {
    // std::cerr << "wrong at sixJ_quick need to set a larger j2a_max" <<
    // std::endl;
    // exit(0);
  } else {
    int k1 = this->JaJbJ_t.vec[j2a][j2b][J2ab];
    int k2 = this->JaJbJ_t.vec[j2b][j2c][J2bc];
    int k3 = this->JaJbJ_t.vec[J2ab][j2c][J2];
    if (k1 < 0 || k2 < 0 || k3 < 0) {
      // cout << "wrong at sixJ_quick k <0 " << endl;
      // cout << "k"
      //      << "\t j2a"
      //      << "\t j2b"
      //      << "\t J2ab" << endl;
      // cout << k1 << "\t" << j2a << "\t" << j2b << "\t" << J2ab << endl;
      // cout << "k"
      //      << "\t j2b"
      //      << "\t j2c"
      //      << "\t J2bc" << endl;
      // cout << k2 << "\t" << j2b << "\t" << j2c << "\t" << J2bc << endl;
      // cout << "k"
      //      << "\t J2ab"
      //      << "\t j2c"
      //      << "\t J2" << endl;
      // cout << k3 << "\t" << J2ab << "\t" << j2c << "\t" << J2 << endl;
      //  exit(0);
    }
    int m = this->TrK_t.exindex2[k1][k2];
    int n = this->TrK_t.exindex3[k1][m][k3];
    int trk = this->TrK_t.vec[k1][m][n];
    sixj = this->vec_6j[trk];
    return sixj;
  }
  return sixj;
}
//
double Pre_coeff::NinJ_quick(int j2a, int j2b, int J2ab, int j2c, int j2d,
                             int J2cd, int J2ac, int J2bd, int J2) {
  //
  if (J2ab < std::abs(j2a - j2b) || J2ab > (j2a + j2b))
    return 0.0;
  if (J2cd < std::abs(j2c - j2d) || J2cd > (j2c + j2d))
    return 0.0;
  if (J2 < std::abs(J2ab - J2cd) || J2 > (J2ab + J2cd))
    return 0.0;
  if (J2ac < std::abs(j2a - j2c) || J2ac > (j2a + j2c))
    return 0.0;
  if (J2bd < std::abs(j2b - j2d) || J2bd > (j2b + j2d))
    return 0.0;
  if (J2 < std::abs(J2ac - J2bd) || J2 > (J2ac + J2bd))
    return 0.0;
  int g2_min = std::max(std::abs(j2b - J2cd), std::abs(j2c - J2bd));
  g2_min = std::max(g2_min, std::abs(J2 - j2a));
  int g2_max = std::min(j2b + J2cd, j2c + J2bd);
  g2_max = std::min(g2_max, J2 + j2a);

  double ninJ = 0.0;
  for (int g2 = g2_min; g2 <= g2_max; g2 = g2 + 2) {
    double part1_t = 0.0;
    part1_t = SixJ_quick(j2a, j2b, J2ab, J2cd, J2, g2);
    if (part1_t == 0.0)
      continue;
    part1_t = part1_t * SixJ_quick(j2c, j2d, J2cd, j2b, g2, J2bd);
    if (part1_t == 0.0)
      continue;
    part1_t = part1_t * SixJ_quick(J2ac, J2bd, J2, g2, j2a, j2c);
    if (part1_t == 0.0)
      continue;
    part1_t = part1_t * phase(g2) * (g2 + 1);
    ninJ += part1_t;
  }
  return ninJ;
}
// ~~~~~~~~~~~~~~~~~~~ Hobra ~~~~~~~~~~~~~~~~~~~ //

double Pre_coeff::G(int &e1, int &l1, int &e_a, int &l_a, int &e_b, int &l_b) {
  // double value_g;
  double part1, part2, part3;
  int l2a = l_a * 2, l2b = l_b * 2, l21 = l1 * 2;
  part1 = std::sqrt((l2a + 1) * (l2b + 1));
  // part2 = this->cg_p.quick(l2a, l2b, l21, 0, 0, 0);

  part2 = this->cg_quick(l2a, l2b, l21, 0, 0, 0);
  //part2 = CGcoeff_2(l2a, 0, l2b, 0, l21, 0);
  // return 1.0;
  if (part2 == 0)
    return 0.0;
  int n2_1, n2_2, n2_3;
  n2_1 = e1 - l1;
  n2_2 = e_a - l_a;
  n2_3 = e_b - l_b;
  // if (n2_1 < 0 || n2_2 < 0 || n2_3 < 0) {
  //   //std::cout << "wrong at HObra G , !! 2n = e - l < 0" << endl;
  //   exit(0);
  // }
  double temp_1 = this->binom_cal(n2_1, n2_2, n2_3);
  double temp_2 =
      this->binom_cal(e1 + l1 + 1, e_a + l_a + 1, e_b + l_b + 1);
  part3 = temp_1 * temp_2;
  // part3 = this->binom_p.vec[n2_1][n2_2][n2_3] *
  //         this->binom_p.vec[e1 + l1 + 1][e_a + l_a + 1][e_b + l_b + 1];
  double value_g = part1 * part2 * std::sqrt(part3);
  return value_g;
  // return value_g;
}

double Pre_coeff::hobra_cal(int &E, int &L, int &e, int &l, int &e1, int &l1,
                            int &e2, int &l2, int &lambda, double &d) {
  double part1;
  double t1, t2;
  if ((E + e) != (e1 + e2))
    return 0.0;
  if (lambda > (l + L) || lambda < std::abs(l - L))
    return 0.0;
  if (lambda > (l1 + l2) || lambda < std::abs(l1 - l2))
    return 0.0;

  int N = (E - L) / 2;
  int n = (e - l) / 2;
  int n1 = (e1 - l1) / 2;
  int n2 = (e2 - l2) / 2;
  int phase_t =
      N + n + n1 + n2; // from the different convension of H.O. Wave_Fun
  int phase_diff = phase(phase_t);

  int e_a, e_b, e_c, e_d, l_a, l_b, l_c, l_d;
  t1 = std::pow(d, (e1 - e));
  t2 = std::pow((1 + d), (e1 + e2));
  part1 = std::sqrt(t1 / (1.0 * t2));
  int e_d_max = std::min(e, e2);
  double part2 = 0.0;
  for (e_d = 0; e_d <= e_d_max; e_d++) {
    e_a = e1 - e + e_d;
    e_b = e - e_d;
    e_c = e2 - e_d;
    double part2_1 = std::pow(-1 * d, e_d);
    double part2_2 = 0;
    for (l_a = e_a; l_a >= 0; l_a = l_a - 2) {
      for (l_b = e_b; l_b >= 0; l_b = l_b - 2) {
        if (l1 > (l_a + l_b) || l1 < std::abs(l_a - l_b))
          continue;
        double G1 = this->G(e1, l1, e_a, l_a, e_b, l_b);
        for (l_c = e_c; l_c >= 0; l_c = l_c - 2) {
          if (L > (l_a + l_c) || L < std::abs(l_a - l_c))
            continue;
          double G3 = this->G(E, L, e_a, l_a, e_c, l_c);
          for (l_d = e_d; l_d >= 0; l_d = l_d - 2) {
            if (l > (l_b + l_d) || l < std::abs(l_b - l_d))
              continue;
            if (l2 > (l_c + l_d) || l2 < std::abs(l_c - l_d))
              continue;

            // double part2_t =
            //     this->ninJ_p.quick(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d *
            //     2,
            //                        l2 * 2, L * 2, l * 2, lambda * 2);
            double part2_t =
                NinJ_quick(l_a * 2, l_b * 2, l1 * 2, l_c * 2, l_d * 2, l2 * 2,
                           L * 2, l * 2, lambda * 2);
            if (part2_t != 0) {
              // double G2 = this->G(e2, l2, e_c, l_c, e_d, l_d);
              // double G4 = this->G(e, l, e_b, l_b, e_d, l_d);
              part2_t = part2_t * G1 * G3 *
                        this->G(e2, l2, e_c, l_c, e_d, l_d) *
                        this->G(e, l, e_b, l_b, e_d, l_d);
            }
            part2_2 += part2_t;
          }
        }
      }
    }
    part2 += part2_1 * part2_2;
  }
  double value = part1 * part2 * phase_diff;
  return value;
}
